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Polarized Superfluidity in the imbalanced attractive Hubbard model 
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We investigate the attractive Hubbard model in infinite spatial dimensions by means of 
dynamical mean-field theory. Using a continuous-time Monte Carlo algorithm in the Nambu 
formalism as an impurity solver, we directly deal with the superfluid phase in the population 
imbalanced system. By calculating the superfluid order parameter, the magnetization, and the 
density of states, we discuss how the polarized superfluid state is realized in the attractive 
Hubbard model at quarter filling. We find that a drastic change in the density of states is 
induced by spin imbalanced populations in the superfluid state. 
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1. Introduction 

The superfluid state in ultracold atomic systems has 
attracted much interest since the successful realization 
of the Bosc-Einstein condensation (BEC) of rubidium 
atoms. 1 In addition to bosonic systems, the superfluid 
state has been observed in two-component fermionic sys- 
tems, 2 where Cooper pairs formed by the attractive in- 
teractions condense at low temperatures. Due to the 
high controllability of the interaction strength and the 
particle number, interesting phenomena have been ob- 
served such as the BCS-BEC crossover 3-5 and the super- 
fluid state with imbalanced populations. 6, 7 These obser- 
vations stimulate further experimental and theoretical 
, investigations on fermionic systems. 

In the existing literature on spin imbalanced popu- 
lations various ordered ground states have been pro- 
posed to be more stable than the polarized superfluid 
(PSF) state, which is naively expected to be realized 
below the critical temperature. One interesting can- 
didate is the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
phase, 8,9 where Cooper pairs are formed with nonzero 
total momentum. This phase has been observed in the 
high field region in CsCoIns, 10-12 and has theoretically 
been discussed in the latter compound, 13 as well as cold 
atoms with imbalanced populations. 14, 15 Another pro- 
posed phase is the breached-pair (BP) phase, where both 
the superfluid order parameter and the magnetization 
are finite at zero temperature. 16-20 When one consid- 
ers higher dimensional optical lattice systems, the BP 
state without momentum dependence may be one of the 
appropriate ground states. It has recently been clarified 
that the PSF state is closely connected to the BP phase 
at half filling in the three-dimensional Hubbard model 
with intermediate attractive interactions. 21 However, the 
Hubbard model has a high symmetry at half filling, 22-24 
and the conclusions may not be applicable to an optical 
lattice system, where the particle density is not fixed at 
half filling due to the existence of the confining potential. 
Therefore, it is important to clarify how the PSF state 
and the BP state are realized in a system away from half 
filling. 



With this purpose in mind, we investigate the attrac- 
tive Hubbard model at quarter filling to discuss the ef- 
fect of the imbalanced spin populations on the super- 
fluid state. By combining dynamical mean-field theory 
(DMFT) 25-28 with the continuous time quantum Monte 
Carlo (CTQMC) method, 29 we study the low tempera- 
ture properties of the system quantitatively. Here, we ex- 
tend the CTQMC method in the continuous-time auxil- 
iary field (CTAUX) formulation 30 to treat the PSF state 
in the Nambu formalism. By calculating the order pa- 
rameter of the superfluid state, the magnetization, and 
the density of states, we clarify the nature of the PSF 
state in the spin imbalanced system. 

The paper is organized as follows. In §2, we intro- 
duce the model Hamiltonian for the attractive Hubbard 
model and briefly summarize the DMFT framework. The 
CTQMC algorithm in the Nambu formalism is explained 
in some detail in §3. In §4, we focus on the attractive 
Hubbard model at quarter filling to discuss how the PSF 
state is realized at low temperatures. A brief summary 
is given in §5. 

2. Model and Method 

We consider a correlated fermion system with attrac- 
tive interactions, which may be described by the Hub- 
bard Hamiltonian, 

where Ci a (cj CT ) is an annihilation (creation) operator of a 
fermion on the iih site with spin a, and n% a = c\ a Ci a . U is 
the onsite attractive interaction, t is the transfer integral 
between sites, [i is the chemical potential, and h is the 
magnetic field. For h = the ground state properties of 
the model have been studied in one dimension, 31-36 two 
dimensions 22,24 ' 37 and infinite dimensions. 23,38-42 Both 
the BCS-BEC crossover and the possibility of a super- 
solid state have been discussed. On the other hand, there 
are few studies addressing the effect of imbalanced popu- 
lations beyond the static mean-field approach except for 
one dimensional system. 15 

To study the infinite dimensional attractive Hubbard 
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model at an arbitrary filling, wc make use of DM FT. 25-28 
In DMFT, the original lattice model is mapped to an 
effective impurity model, which accurately takes into 
account local particle correlations. The lattice Green's 
function is obtained via a self-consistency condition 
imposed on the impurity problem. This treatment is 
formally exact in infinite dimensions, and the DMFT 
method has successfully been applied to strongly cor- 
related fcrmion systems. 

When the supcrfluid state is treated in the framework 
of DMFT, the local self-energy should be described by a 
2x2 matrix as 



f /. \ f S t (iw n ) S(iu) n ) 



(2) 



where Yi a {iw n ) [S^iu^)] is the diagonal (off-diagonal) el- 
ement of the self-energy in the Nambu formalism and 
the Matsubara frequency is uj n = (2n + l)ir/f3, with (3 
the inverse temperature. Note that we do not take into 
account /c-dependent correlations, but dynamical correla- 
tions through the frequency-dependent self-energy. This 
enables us to discuss the stability of the s-wave superfluid 
state more quantitatively beyond the static mean-field 
theory. 

The lattice Green's function is then given in terms of 
the local self-energy as, 

G~ 1 {k,iui n ) = (iuj n + h)a + (fi- e k ) cr z ~t> (iui n ) , (3) 

where ctq and ct z are the identity matrix and the z- 
component of the Pauli matrix, and is the dispersion 
relation for the non-interacting system. The local lattice 
Green's function is obtained as, 



G{ioo n ) = J dkG(k 7 iui n ). 



(4) 



In the calculations, we use the semi-circular density of 
states, p(x) = l/Nj2k S ( x ~ £ k) = 2/ttD^I - {x/D) 2 , 
where D is the half bandwidth. The self-consistency 
equation 43 is then given by 



D 



a z G(iu!„)a z 



(5) 

When one discusses low energy properties in strongly 
correlated systems in the framework of DMFT, an impu- 
rity solver is necessary to obtain the Green's function and 
the self-energy for the effective impurity model. There 
are various numerical techniques such as exact diagonal- 
ization 44 and the numerical renormalization group. 45-47 
A recently developed and particularly powerful method 
is CTQMC. In this method, Monte Carlo samplings of 
collections of diagrams for the partition function are per- 
formed in continuous time, and thereby the Trotter er- 
ror, which originates from the Suzuki- Trotter decompo- 
sition, is avoided. Furthermore, this method is applica- 
ble to more general classes of models than the Hirsch- 
Fye algorithm. 48 The CTQMC method has successfully 
been applied to various systems such as the Hubbard 
model, 49 ' 50 the periodic Anderson model, 51 the Kondo 
lattice model 52 and the Holstein-Hubbard model. 53 



3. Continuous-Time Quantum Monte Carlo sim- 
ulations in the Nambu Formalism 

In this section, we explain the CTAUX method, 30 and 
extend it to treat the supcrfluid state. A similar solver 
was recently proposed, 51 where the supcrfluid state is 
treated by means of a canonical transformation. The An- 
derson impurity model we have to solve is given by 

H = Ho + Hu, (6) 
H = ^2 ep<jn pa + ^ (Vpad^apa- + h.c.) 



pa 



pCJ 



E ( A AA 



pi 



h.c 



H v = -U 



rid^ndi - - (n d f + ridi - 1) 



^E dcr n d<T , (7) 
(8) 



where a pa (d a ) annihilates a fermion with spin a in the 
pth orbital of the effective baths (the impurity site). 
e p(T and A p represent the effective bath, and V pa rep- 
resents the hybridization between the effective bath and 
the impurity site. E da is the energy level for the impu- 



rity site, n pa = a vrT a pa , 



and n da = d\d a . We note that 



the total particle number is not conserved in the model. 
The Green's functions should be defined by G(t) = 
(T T ^(T)^ t (0)), where T T is the imaginary-time ordering 
operator and ipH T ) = ( C \( T ) c i( T ))- The Green's func- 
tions are 2x2 matrices with elements 



G(t) 



G t (r) F(r) 
F*(t) -G;(-t) 



(9) 



where G a (r) — (T T c cr (r)cJ.(0)) denotes the normal 
Green's function, and F(t) = (T t c^(t)cj_(0)) and 
F*( T ) — (^Tcj(r)cj(O)) anomalous Green's functions. 
Here, we have chosen the Green's functions G a {r) to be 
positive. 

To perform simulations, we consider here a weak cou- 
pling CTQMC approach. The partition function Z is 
given by 

Z = Tr\e-^T T e-fo d rH 2 (r) 

OO 

- E 

n=0 

x Tr[e- (/3 - r " )Hl (-#2)e 



dn / dr 2 - ■ ■ dr n 

Jti Jl„_i 



(t„-t„_i)Hi 



- { - T2 - Tl)Hl {-H 2 )e 



-T 1 H 1 



(10) 



where we have divided the impurity Hamiltonian Eq. (6) 
into two parts as, 



Hi 
H 2 



H-H 2 , 

Hu - K/P 
K 



g7s(«T+"4-- 1 ) j 



(11) 



(12) 



with 7 = cosh _1 (l + (3U/2K), and K some nonzero 
constant. The introduction of the Ising variable 
s in H 2 enables us to perform simulations away 
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from half-filling. An nth order configuration c = 
{sii S2, ■ ■ ■ , s n ; Ti, T2, ■ • • , r n } corresponding to auxiliary 
spins Si, S2, . . . , s n at imaginary times t\ < T2 < . . . < t„ 
contributes a weight 



(13) 



to the partition function. Here, Z = Tr[e- pHl ] and 

is an n x n matrix, where each element consists of a 2 x 2 

matrix: 



f(n) 
ij 

' ij 

An) 
3ij 



p(") _ g( n ) (™) _ /(™)^ ; 



5ot( r i 



t;) /o(t; 



-n) 



(14) 

(15) 
(16) 

,(17) 



with i, j = 1, 2, • • • n. 

The sampling process must satisfy ergodicity and (as a 
sufficient condition) detailed balance. For ergodicity, it is 
enough to insert or remove the Ising variables with ran- 
dom orientations at random times to generate all possible 
configurations. To satisfy the detailed balance condition, 
we decompose the transition probability as 



(18) 



where pP ro P(p acc ) is the probability to propose (accept) 
the transition from the configuration i to the configura- 
tion j. Here, we consider the insertion and removal of the 
Ising spins as one step of the simulation process, which 
corresponds to a change of ±1 in the perturbation order. 
The probability of insertion/removal of an Ising spin is 
then given by 



P 



prop 



{n — > n + 1) = — 



dr 

2/3' 



P 



prop 



(n + 1 — > n) 



1 



1 



(19) 
(20) 



For this choice, the ratio of the acceptance probabilities 
becomes 

dct JVW 



pace ^ n _^ n _|_ -Q 



A" 



-7S»+1 . 



1 



det AT("+!) 



(21) 



pace + 1 -> n) 

When the Metropolis algorithm is used to sample the 
configurations, we accept the transition from n to n ± 1 
with the probability 



nil) 



(22) 



' pace (ti ± 1 — ^ n ) 

In each Monte Carlo step, we measure the following 
Green's functions (0 < r < (3), 



G„(t) = 
F(r) = 
F*(t) = 



iTr[e-^ C(7 (r)ct(0)] 
iTr[ e -^c t (r)c;(0)] 



iTr 

z 



e-^ci(r)4(0) 



(23) 
(24) 
(25) 



configuration c is given by 
G c a (T) = det[A( n) ]det 



F c (t) = det[AfW]det 
F* c (t) = det det ( [ N< £\ 



r/v(")l _1 


n 


JX<j 


SOo-(t) 




Q' 


R' 


fo(r) 




Q*' 


R*' 


ft(r) 



(26) 
(27) 
(28) 



where Q a ,Q' ,Q*' , R&, R' , Ft*' are vectors, in which the 
ith clement (i = 1, 2, • • ■ , n) is defined by 

Qti = {-SOl-fa) /o( T 0} T , 



<5l* = {/o(t - r) 5o4.(t - ^)} T , 



By using Wick's theorem, the contribution of a certain 



(29) 
(30) 

Qi = i-Mn) -goi(-n)} T , (31) 
Q*' = Qti. (32) 
Rn = (e 7Sl - l){flbt(r - t0 /o(r -T*)}, (33) 
% = (e 7 * - l){/o(-ri) -.904(^)1, (34) 

i?- = (35) 

= ( e 7^_l ){/ *( T _ Ti) _^ (r ._ T ) } .( 36 ) 

In this paper, we use the half bandwidth D as the unit 
of the energy and set K = 1 in the CTQMC simula- 
tions. We thus calculate static physical quantities such 
as the order parameter of the supcrfluid state A and the 
magnetization m, which are defined by 

A = (c t c L )=F(0 + ), (37) 
m = 5>(4O = -5> G '( +)- (38) 

a a 

Furthermore, by applying the maximum entropy method 
(MEM) to the Green's functions, we deduce the spectral 
functions, which allows us to discuss static and dynami- 
cal properties of the system. 

In Fig. 1, we show, as an example, the normal and 
anomalous Green's functions when U = 1, h = 0.1 and 
T = 0.01. The Green's functions were measured on a 
grid of a thousand points. In this case, the system has 
both a magnetization m ~ 0.005 and a superfluid order 
parameter A ~ 0.2. Therefore, we can say that the PSF 
state is realized in this parameter region. Note that a 
large difference appears between Gf(r) and Gj.(t) near 
t ~ and P although the magnetization is small. This 
may affect dynamical properties. 

4. Superfluid state in a magnetic field 

Here, we focus on the attractive Hubbard model at 
quarter filling to discuss how the PSF state is realized 
at low temperatures. First, we perform calculations at a 
fixed temperature. Results for the systems with weak (in- 
termediate) coupling [U — 1 (U = 2)] are shown in Fig. 
2. When no magnetic field is applied, the system is in 
the supcrfluid state at low temperatures. In fact, we find 
that the supcrfluid gap opens around the Fermi level and 
that peak structures appear at the edges of the gap in 
the density of states, as shown in Fig. 3. These results are 
consistent with those obtained by other groups. 40 ' 42 If a 
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Fig. 1. Green's functions as a function of r//3 in the quarter- 
filled system at U = 1, h = 0.1 and T = 0.01. The solid (dashed) 
line represents the Green's function for the up (down) spin and 
the thin line the anomalous Green's function. The inset shows 
the probability distribution for configurations with perturbation 
order n at the temperatures T = 0.01,0.014 and 0.025. 



magnetic field is applied to the system, these peaks move 
to low (high) energy region in the density of states for 
up (down) spin. Pairing correlations are then suppressed, 
and a magnetization is induced, as shown in Fig. 2. We 
note that at low temperatures, the introduction of a mag- 
netic field has little effect on the static quantities A and 
m, but produces a drastic change in the density of states. 
In fact, it is found that when U = 2 and h = 0.3, one of 
the peaks disappears and the other remains above (be- 




Fig. 3. Solid (dashed) lines represent the spectral functions for 
fermions with up (down) spin when U = 1, T = 0.04 (a), and 
U = 2,T = 0.05 (b). 




Fig. 2. The superfiuid order parameter and the magnetization as 
a function of the magnetic field when U = 1,T = 0.04 (a), and 
U = 2,T = 0.05 (b). The insets show the critical behavior for 
the order parameter. 



low) the Fermi level in the density of states for up (down) 
spin although the superfiuid gap is still open. Therefore, 
we conclude that dynamical properties are strongly af- 
fected by the spin unbalanced populations. A further in- 
crease in the magnetic field smears the superfiuid gap 
around the Fermi level and the superfiuid order parame- 
ter vanishes. This suggests the existence of a phase tran- 
sition to the normal metallic phase. By examining the 
critical behavior A ~ \h— h c \ x / s with the exponent (5 = 3, 
we obtain the critical fields h c (U = 1, T = 0.04) - 0.0885 
and h c (U = 2,T = 0.05) ~ 0.425, as shown in the in- 
sets of Fig. 2. It is also found that the phase transition 
induces a cusp singularity in the magnetization curve. 
The results obtained here are in contrast to those in the 
half-filled attractive Hubbard model on the simple cubic 
lattice, where the PSF state smoothly connects to the 
normal metallic phase. 21 This may result from the fact 
that the competition between the superfiuid state and 
the charge density wave state enhances fluctuations for 
the superfiuid order parameter due to the high symme- 
try at half filling. It would be interesting to clarify this 
point, which is beyond the scope of our study. 

We also show the temperature dependence of the su- 
perfiuid order parameter in Fig. 4. When h = 0, as tem- 
perature is decreased, the order parameter A appears 
where the phase transition occurs from the normal metal- 
lic state to the superfiuid state. By examining the critical 




behavior A ~ \T — T C \P with the exponent j3 = 1/2, we 
obtain the critical temperature T c ~ 0.095, as shown in 
the inset of Fig. 4. On the other hand, when the magnetic 
field is switched on, pairing correlations are suppressed. 
For h = 0.3, it is found that the superfiuid order param- 
eter is decreased and the critical temperature is shifted 
to T c ~ 0.077. A large magnetic field destroys the super- 
fluidity and the normal metallic state is realized instead. 
In fact, we could not find any finite A down to low tem- 
peratures (T = 0.033) in the case h = 0.6. Note that the 
two curves saturate almost at the same value of A when 
T — > 0. This suggests that by increasing the magnetic 
field at zero temperature, the superfiuid ground state is 
little affected and eventually a first order phase transition 
occurs to the normal metallic state. Therefore, it may be 
difficult to realize the BP state with finite magnetization. 

To clarify this, we next examine how the magnetization 
appears in the superfiuid state. In Fig. 5, we show a semi- 
log plot of the magnetization normalized by the applied 
field. When a tiny magnetic field is applied to the system, 
m/h corresponds to the magnetic susceptibility %. In the 
nonintcracting system (U = 0), m/h saturates at low 
temperatures at the value x(T = 0) = 2p(x = 0) ~ 1.16, 
and in the interacting case, our results are consistent 
with those obtained by Keller et al. 39 Increasing the 
attractive interaction in the presence of a finite mag- 
netic field, fermion pairs are formed at high temperatures 



m/h 




and thereby magnetic correlations are suppressed and 
the magnetization decreases. When the magnetic field is 
small enough, a phase transition occurs to the superfiuid 
state at low temperatures. In this state, the magnetiza- 
tion rapidly decreases below the critical temperature, as 
shown in Fig. 5. This means that it is difficult to real- 
ize the BP ground state with finite A and m at zero 
temperature. To confirm this, we also show the quantity 
At(= — Tlogm) as a function of temperature in Fig. 6. 
When T — > 0, the data approach a finite value in the 
superfiuid state, while they approach zero in the normal 
metallic state. This means that the magnetization decays 
exponentially in 1/T in the superfiuid state. Therefore, 
we can say that the BP state is not realized in the ground 
state, at least, in this quarter-filled system. 

By performing similar calculations, we have obtained 
the phase diagram for the spin imbalance parameter 
P[= (iV t — N±) j (iV-f + Ni) = 2m], which is sometimes 
used in the discussion of optical lattice systems, as shown 
in Fig. 7. When the temperature decreases with fixed im- 
balanced populations, a phase transition occurs to the 
PSF state. Figure 8 shows the density of states for each 
spin component in a system with P ~ 0.02 and U = 2. 
It is found that at high temperatures (T > T c ), the nor- 
mal metallic state is realized, where the spin imbalanced 



0.05 




Fig. 5. Normalized magnetization m/h as a function of temper- 
ature T in the system with U = 0,1 and 2. 



Fig. 7. Phase diagrams for the quarter-filled system with (7 = 1 
and U = 2. Open (solid) circles indicate the normal (PSF) state 
and phase boundaries are guides to eyes. 
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Fig. 8. Density of states for the system with {7 = 2 and P = 0.02 
when T = 0.2, 0.1, 0.05 and 0.033. 



populations have little effect on the density of states. By 
contrast, in the superfluid state (T < T c ), a large dif- 
ference appears in the low energy region of the spectral 
functions. Since the populations with up spin are slightly 
larger than those with down spin in the case considered, a 
certain energy is necessary to add a fermion with up spin 
in the superfluid state. Therefore the low energy peak ap- 
pears above the Fermi level in the density of states for 
up spin at low temperatures. As temperature is lowered 
to zero, it may be difficult to realize a state with P ^ 0, 
as discussed before. We cannot rule out that a phase 
transition from the superfluid phase back to the metal- 
lic phase (reentrant behavior) will occur at temperatures 
below the range accessible to us. In any event, the imbal- 
anced populations should play a crucial role at very low 
temperatures, in particular, in the dynamical properties. 

5. Summary 

We have investigated the attractive Hubbard model in 
infinite dimensions by means of DMFT. Here, we have 
used the CTAUX method as an impurity solver, which 
has been extended to treat the superfluid state directly in 
the Nambu formalism. We have calculated the superfluid 
order parameter, the magnetization, and the density of 
states systematically to discuss how the PSF state is re- 
alized at low temperatures. It was found that when the 
temperature is lowered in the presence of a fixed mag- 
netic field, a superfluid phase transition indeed occurs in 
our model, and the magnetization exponentially decays 
in the superfluid state. This suggests that the BP phase is 
unstable at zero temperature. We have also found that a 
drastic change in the density of states is induced by spin 
imbalanced populations in the superfluid state although 
the spin imbalance has little effect on static quantities. 
It is an interesting problem to clarify how such dynam- 
ical properties are realized in a low dimensional optical 
lattice with a confining potential, which is now under 
consideration. 
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